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Abstract 


A new  quadrilaterial  finite  element  for  the  problem  of  bending  of  a 
uniform,  homogeneous,  isotropic,  flat  plate  having  an  arbitrary  number  of 
nodal  points  and  allowing  an  arbitrary  level  of  refinement  is  developed 
from  solution  functions  to  the  governing  differential  equation.  Interelement 
continuity  to  the  appropriate  level  of  normal  derivatives  is  provided  in 
least  error  fashion  by  the  use  of  spline  functions  for  the  edge  displacement 
quantities  and  minimization  by  the  Galerkin  technique.  Even  for  general 
quadrilateral  domains  no  numerical  quadrature  is  required  for  developing 
this  finite  element  since  all  required  integrals  are  taken  over  straight 
boundary  edges  and  can  be  evaluated  explicitly.  With  this  new  type  of  element, 
an  improvement  in  accuracy  is  expected  to  result  when  large  quadrilateral  regions 
are  modeled  by  a single  element  with  an  appropriate  number  of  edge  nodal 
points,  because  the  error  of  this  finite  element  approximation  is  associated 
with  the  boundary  region  of  the  elements.  The  force  vector  corresponding 
to  the  transverse  pressure  is  derived  from  a series  solution  to  the  inhomo- 
geneous differential  equation,  so  that  the  effect  of  the  forcing  function 
in  the  element  interior  can  be  evaluated  to  any  degree  of  accuracy  desired. 
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1,  Introduction 


The  vast  majority  of  boundary  value  problems  do  not  lend  themselves 
to  solution  by  an  exact  method.  Analysts  therefore  resort  to  approximate 
methods  that  reduce  the  continuum  problem  to  the  problem  of  finding  a finite 
number  of  parameters  that  characterize  the  solution.  The  most  versatile 
approximation  method  to  date  for  solving  linear  boundary  value  problems  is 
the  finite  element  method  (FBi).  This  method  is  fully  automatible  and  pro- 
duces in  general  sparcely  populated  matrices  which  lend  themselves  to  inversion 
by  techniques  that  are  both  economic  in  storage  and  saving  in  computer  time, 
thus  permitting  application  to  large  systems.  Although  the  FBI  has  initially 
been  developed  by  structural  analysts  and  most  of  its  vocabulary  is  derived 
from  structural  mechanics,  its  versatility  has  been  recognized  and  its  appli- 
cation has  spread  to  other  fields  of  continuum  mechanics. 

The  FBI  is  still  a subject  of  vigorous  research  by  both  engineers  and 
mathematicians  alike.  Efforts  are  going  on  to  enlarge  the  field  of  applica- 
tion of  the  FBI,  to  give  information  about  the  error  in  the  solution,  and 
to  improve  upon  existing  finite  element  approximations  to  continuum  problems. 
These  efforts  include  reformulating  the  coefficient  matrix,  which  in  finite 
element  nomenclature  is  called  stiffness  matrix,  in  such  a manner  as  to 
achieve  reduced  computer  time,  refined  representation  of  the  boundary,  and 
improved  accuracy.  This  work  is  part  of  that  effort.  The  point  of  depart- 
ure for  this  present  work  is  the  extended  field  method  (EXFM) , a method 
of  analysis  that  exhibits  remarkable  convergence  but  lacks  the  flexibility 
of  application  of  the  FEM.  In  this  report  (dissertation),  the  EXFM  is  used 
to  formulate  a new  finite  element  stiffness  matrix  and  the  associated  force 
vector  for  the  problem  of  plate  bending. 
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In  brief,  this  new  formulation  employs  solution  functions  in  the  interior  of 
the  finite  element  and  finds  some  best  fit  along  the  element  boundary. 

Inter-element  compatibility  in  the  sense  of  P ian  (ref.  1 .)  is  achieved 
by  using  a Galerkin-type  error  reduction  along  the  boundary.  Instead  of 
deriving  the  force  vector  that  corresponds  to  the  forcing  function  acting 
over  the  interior  of  the  element  from  work  equivalence,  as  is  customary 
in  standard  finite  element  foimulations,  the  force  vector  is  here  derived 
from  exact  solutions  to  the  inhomogeneous  differential  equation.  The 
stiffness  kernel  is  obtained  from  virtual  work  considerations  that  require 
integration  of  the  stress  resultants  along  the  element  boundary.  There 
is  no  need  for  integration  over  the  domain  of  the  finite  element,  since 
all  displacement  functions  employed  in  this  formulation  satisfy  the 
differential  equation  exactly.  This  amounts  to  an  equilibrium  formu- 
lation. Choosing  for  the  solution  functions  certain  series  whose 
recurrence  relation  is  simply  determined  by  an  integer  that  multiplies  the 
argument,  this  new  finite  element  can  be  easily  refined  by  increasing  both 
the  number  of  displacement  functions  and  the  number  of  nodal  points  along 
the  element  boundary.  This  process  can  be  automated.  Tliis  flexibility  in 
the  number  of  nodal  points  per  element  permits  the  analyst  to  use  a suffi- 
ciently refined  element  of  this  type  instead  of  an  equivalent  substructure 
consisting  of  a number  of  finite  elements.  Instead  of  progressive  parti oning 
the  researcher  can  use  the  process  of  progressive  element  refinement 
for  demonstrating  convergence.  It  is  noted,  that  for  an  equivalent  refinement  of  the 
customary  finite  element  solution  the  growth  of  degrees  of  freedom  of  the  system 
is  roughly  proportional  to  the  square  of  the  number  of  edge  partitions  per 
element  edge  of  the  original  element  size  in  the  case  of  progressive  par- 
titioning, while  for  element  refinement,  as  is  possible  for  this  new  formu- 
lation, the  number  of  degrees  of  freedom  grow's  proportional  to  the  partitioning 
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of  the  circunference.  (Starting  with  a four  node  element,  the  proportionality  fac- 
tors are  i)  for  progressive  partitioning:  S,  “4  , 9 , 5^  25 

^ and  ii)  for  element  refinement:  5,  = 2!*  ), 

2.  The  Plate  Bending  Finite  Element  Derived  by  the  Direct  Technique 
2.1  The  direct  technique  for  the  plate  bending  finite  element 

IVhile  the  techniques  now  generally  employed  for  ihe  formulation  of  finite 
element  stiffness  matrices  are  based  upon  the  method  of  variational  calculus  (Ref's 
2,3,4),  the  approach  taken  here  is  the  direct  technique,  or  as  it  is  tenrtdin  some 
literature  the  equilibrium  technique.  The  general  approach  of  the  equilibrium  tech- 
nique is  to  obtain  exact  solutions  to  the  governing  differential  equation  that  are 
expressed  in  terms  of  displacement  quantities  at  the  nodal  points,  equate  the  boun- 
dary nodal  point  displacement  quantities  of  the  adjacent  elements  and  enforce  equili- 
brium of  the  force  quantities  at  the  nodal  points.  For  one- dimensional  continua  such 
as  rods  and  beams  this  procedure  is  simple  to  apply,  as  the  total  boundary  consists 
of  two  discrete  points  and  the  stress  resultants  are  already  discretized  quantities. 
In  the  case  of  two  dimensional  problems,  such  as  the  plate  bending  problem,  proper 
continuity  of  the  primary  variable  across  the  entire  element  boundary  cannot  always 
be  achieved  and  equilibrium  can  only  be  achieved  in  some  average  sense.  The  method 
described  in  detail  later  in  this  report  allows  discontinuities  of  the  primary  vari- 
able across  element  boundaries  but  minimizes  these  discontinuities  in  the  Galerkin 
sense,  and  enforces  equilibrium  in  an  average  sense  by  appeal  to  the  priniciple  of 
virtual  work. 

To  develop  this  approach  for  the  case  of  a plate  bending  finite  element,  con- 
sider the  transverse  displacementW(#][rhe  transverse  displacement  satisfies 

D v'^w  - p(x) , 

where  pCX)  , the  forcing  function,  is  the  transverse  pressure  and  u is  the  bending 
rigidity  of  the  uniform  plate.  The  transverse  displacement  may  be  expressed 
as  the  sum  of  a solution  to  the  inhomogeneous  equation,  Wp(X),  satisfying 
come  convenient  boundary  condition  and  a solution  of  the  homogeneous  equation^ 


w,(x) , 


with  adjustable  coefficients  that  make  it  possible  to  achieve 
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any  boundary  displacement  within  the  bounds  of  some  required  accuracy.  Summing 
these  components  the  transverse  displacement  is  given  by 


W(X)  - Wp(X)  + W^(X)  (2) 

The  solution  to  the  homogeneous  equation  includes,  of  course,  both  elastic 
as  well  as  rigid  body  displacements. 

The  development  of  the  equilibrium  finite  element  requires  a finite 
number  of  displacement  functions  that  are  dependent  upon  undetermined  para- 
t meters.  Truncated  series  solutions  are  selected  for  the  displacement 

functions,  and  they  may  be  written  in  the  form 

N M 

! Wfc(X)  = I and  Wp(X)=I  B(X)p^  , 

<-l  4- 1 

where  the  p.j  are  Known  pressure  parameters  that  are  derived  from  the  forcing 
function  and  the  h-j  are  the  undetermined  parameters.  The  displacement 

^ • 

boundary  conditions  for  the  plate  bending  problem  are  the  prescription  of 
w(£)  , the  displacement,  and  of  W,n(S)  the  rotation  normal  to  the  boundary, 
where  is  a boundary  coordinate.  The  corresponding  nodal  degrees  of  free- 

dom at  a given  point  onthe  boundary  are 

qj  = 

qj..  * 

Here  the  comma  behind  the  variable  followed  by  an  index  quantity  that  is  a 
spatial  variable  represents  partial  differentiation  of  the  dependent  variable 
with  respect  to  that  spatial  variable  and  n is  the  local  normal  at  the  boundary. 

A transformation  from  the  nodal  degrees  of  freedom,  the  q^tto  the  con- 
tinuous parameters,  the  h^,must  be  obtained  so  that  the  displacement  functions 
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can  be  expressed  in  terms  of  nodal  quantities.  One  possible  approach  is 
to  obtain  the  inverse  transformation  by  simply  evaluating  the  nodal  dis- 
placement quantities  in  terms  of  the  displacement  functions: 

q, " w {«,)  ■ t h,  + i 

q,  - - I + I 

q,-  w,,(e,)-  i,  Jf4,,(4.)h,  + [ 

where  k - N/3  . 

In  matrix  form  the  above  can  be  written  as 

q “ B hi  +-  C p . (3) 

For  the  purpose  of  obtaining  the  in  terms  of  it  is  required  that  N“N 
so  that  B is  a square  matrix.  This  transformation  enforces  continuity 
of  the  displacement  quantities  at  the  nodal  points  but  says  nothing  about 
inter-element  continuity  of  the  displacement  quantities  along  the  element 
boundary  away  from  the  nodal  points.  For  this  work  a more  complex  approach 
is  taken  that  improves  upon  the  match  of  the  displacement  quantities  of 
adjacent  elements  along  their  common  edges.  The  edge  displacement  quanti- 
ties W(^)  and  are  expressed  directly  in  terms  of  the  nodal  degrees 

of  freedom  * 

N 

w(4)  * Ww  ■ Z 5i(4)q; 

vy?)  -f  (4)  - Z R.(4)q. 

in  such  a manner  that  the  nodal  displacement  quantities  on  a comnron  edge  of 
two  adjacent  elements  determine  the  same  displacement  and  rotation  of  this 
edge  for  both  elements.  The  discrepancy  between  this  exterior  representation 
of  the  displacements  and  rotation  and  the  representation  of  these  quantities 
by  the  solution  functions  for  the  element  interior  is  then  minimized. 
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Thus,  interelement  compatibility  in  the  sense  of  Pian  (ref.l)  is 
achieved.  Perhaps  the  most  useful  edge  displacement  functions  W(^) 
and  are  cubic  splines  for  the  transverse  displacements  of  the 
edge  and  linear  splines  for  the  edge  rotation.  Such  an  element  would 
be  interelement  compatible  with  some  of  the  finite  elements  currently 
in  use,  such  as  the  QUAD4  element  of  NASTRAN  or  the  XASTRAI4  beam  element, 
which  have  this  displacement  pattern.  These  boundary  displacement 
functions  would  also  allow  assemblies  of  finite  elements  \diere  a comer 
nodal  point  of  one  element  coincides  with  a mid-edge  nodal  point  of 
another  without  violating  the  interelement  compatibility  requirement. 

The  transformation  from  the  nodal  degrees  of  freedom,  the  , to  the 
distributed  parameters,  the  hj  . is  then  obtained  by  requiring  that  the 
discrepancies  of  the  displacement  quantities  between  the  exterior  and  the 
interior  representation  be  orthogonal  to  some  set  of  independent  test 
functions.  With  the being  such  a set  of  test  functions  and  the  dis- 
crepancies given  by 

error  Jw(?)]  = w(?)-w(?) 
error  (w,.(«)]  - - f (?)  , 

the  transformation  obtained  from  the  orthogonality  requirement 

f error  =0  , 


is  given  by  the  matrix  equation 

Adj'Bhi  -t-Cp 

The  upper  index  limits  must  be  chosen  so  that  the  coefficient  matrix 


f4) 
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of  the  distributed  parameters  will  be  square.  The  inversion  proper- 
ties of  the  coefficient  matrix  B depend  strongly  upon  the  choice  of 
the  set  of  test  functions.  In  particular  none  of  the  test  functions 
may  be  orthogonal  to  the  set  of  component  functions  of  the  error.  Similarly 
none  of  the  component  functions  of  the  error  may  be  orthogonal  to  the 
set  of  test  functions.  If{(P^(^)}  is  a set  of  test  functions  satisfying 
these  criteria,  then  the  inverse  of  the  coefficient  matrix  B exists. 

To  minimize  the  error  that  results  from  numerical  inversion  it  is 
furthermore  desirable  that  the  [(^j be  an  orthogonal  set  rather  than 
merely  being  linearly  independent.  The  choice  of  this  set  of  ortho- 
gonal test  functions  is  discussed  further  in  sections  2.3  and  2.4a. 

Having  obtained  the  transformation  from  nodal  degrees  of  freedom  to  the 

parameters  of  the  solution  function  by  multiplying  equation  (4)  with  ft 

h = - B’'£p  , 

it  is  now  required  to  obtain  the  stiffness  kernel.  In  a variational 
approach  it  would  be  necessary  to  integrate  the  strain  energy  inte- 
gral over  the  two-dimensional  domain.  For  the  general  quadrilateral 
element  this  would  require  numerical  quadrature.  In  the  equili- 
brium approach  the  displacement  functions  chosen  already  satisfy 
the  differential  equation,  therefore  there  is  no  need  to  apply 
the  Principel  of  the  Stationary  Total  Potential  . Instead  it  is 
required  to  enforce  equilibrium  at  the  boundary.  An  error  orthogonal i- 
zation  scheme  similar  to  the  one  used  for  the  displacement  functions, 
where  the  discrete  force  quantities  at  the  nodal  points  are  treated  in 
the  distributional  sense  will  in  general  not  work  if  the  set  of  test 
functions  is  an  orthogonal  set.  This  becomes  evident  if  one  considers 
an  element  where  all  nodal  points  are  located  at  the  zeros  of  one  of 
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the  test  functions.  The  corresponding  integral  would  yield  a row  of  zeros 
for  the  coefficient  matrix  and  render  it  singular.  Using  a set  of  merely 
independent  test  funct*  ns  that  would  circumvent  this  problem  is  deemed 
not  satisfactory  when  primarily  accuracy  is  sought. 

A rational  approach  is  to  employ  the  Principle  of  Virtual  Work  to 
obtain  a relation  between  the  generalized  nodal  forces  and  the  stress 
resultants  al.ong  the  boundary  of  the  finite  element.  Briefly  stated,  write 

m.  - /M, . 

Here  (ref. 5) 

^ “ 4*  i/n  ~ ^ Mr)  6%  + Z(Mn6“Mr)i)<iW|i 

k«l 

is  the  virtual  work  done  by  the  stress  resultants  and  comer  forces 
while  going  through  the  boundary  edge  virtual  displacement  and 


is  the  virtual  wcrk  done  by  the  generalized  nodal  forces  while  going 
through  the  nodal  virtual  displacements.  The  quantities  Mn,Mrisand  Vi. 
in  equation  (7  ),  i.e.  the  normal  moment,  the  twisting  moment  and  the 
Kirchhoff  shear  respectively  at  the  plate  edge  are  linear  combinations 
of  partial  derivatives  with  respect  to  the  local  edge  coordinates  (ref  5) 
Using  appropriate  transformation  coefficients  these  quantities  can 
be  expressed  in  terms  of  the  reference  coordinates  (x,y). 

Thus,  from  Fig.  (5) 

V„  - - D 
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which  when  substituted  into  the  virtual  work  expression  give  for 
the  polygonal  plate,  (Fig.  1 ) 


Z {t„-  [W,„(>U.)-W,„(X„j] 

+ tin'  W,»,(X«J  <Sw(X»J 
>•» 

+ Jw<?)  d% 

+ i„  f .!w(5)  d? 
tlskf  w,.^(e)<Sw(?)  de 
k <«■  f w.toC?)  (!ww  d? 
+ i,,' ^ W,  **(^)  dw,x{^)d^ 
+ le,'  I W,5(^(£) 


+ !«■)  w,^(?l  dw, ,(?)<!« 

+ -tok  w,«(?)  (5w,,(5)  dft 

+ tik'  S.  W,,j(5)  (iw.j^)  ds 
+ t«r  iw„(«)d^  _■ 
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(10) 
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where 

-D(i-D)[«t-(iC+(^U 

i-  -D  l(^i+  (2-p)(^)ii):  ] 


-D  [(l-2u)(^a-(2- 

ts" 

"D  [ (l-2D)(^)ii(^)k-(2- 

-D  i($);  + (2-u)(^i^):i 

0 l&,  + 1 

isr 

-D 

D l&A  + 

-D 

ink 

D 2(1 -P)  «), 

Substituting  for  W and(jw,  the  virtual  work  expression  may  be  cast 
in  matrix  form 

(SfF  = ill)  (fill)  +Rpi}  «« 

Substituting  the  transformation  equation,  eq.  (5)  into  the  virtual  work 
expression,  eq.  (11)  gives 

lifF  - -A'B''((DBX-R)p} 


II 


Quadrilateral  Plate 
Fig.l 


which  holds  for  arbitrary 


so  that  the  equilibrium  equation 


takes  the  form 


where 


F ' IK^i  -F  , 
IK 


and 


The  required  coefficient  matrices  in  the  case  of  a quadrilateral  plate 
bending  finite  element  are  developed  in  the  sequel.  The  components  of 
the  matrices  /A  ^ S (C  are  the  corresponding  coefficients  in  equa- 
tions (22a)  and  (2Cb)  section  2.4a  and  the  vectors  h) , (I|  and  p are 
defined  in  section  2.4b  below.  The  components  of  the  vector  F corres- 
pond to  those  of  the  vector  (I|  in  the  sense  that  their  inner  product 
is  work.  As  a result  of  substituting  the  series  representations  of  the 
displacement  function  into  equation  (10)  the  matrices  (D  and  R are  each 
linear  combinations  of  matrices  defined  by  equation  (23)  in  section  2.4b 
below. 

2 . 2 The  Displacement  Functions  for  the  Plate  Bending  Finite  Element 
An  obvious  choice  for  the  solution  function  Wf(%)  for  the 
inhomogeneous  equation  is  provided  by  the  Navier- solution  vdiich 
satisfies  simply  supported  boundary  conditions  at  the  edges  of  a 
rectangular  reference  field.  For  the  homogeneous  biharmonic  equation 
the  complementary  solution  could  be  chosen  to  be  a polynomial 
consisting  of  a complete  set  of  linearly  independent  polynomial 
terms  up  to  a given  order  that  satisfy  the  homogeneous  equation 
identically.  Such  a set  of  solution  functions  up  to  the  6th  order  is 


x’ 

A'' 

^ M 

^ l ^ r r 

% i t s r„v.»  ic-.*!. 
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It  is  not  difficult  to  show  that  for  the  subset  of  purely 
order  polynomial  terms  there  exist  at  most  four  linearly  independent 
solutions  to  VvV“0  and  if  n ^ 3 there  exist  exactly  four. 

Considering  the  simpler  transformation  equation,  eq.  (3  ),  end  a 
polynomial  series  of  solution  functions  to  the  homogeneous  biharmonic 
equation,  one  will  find  it  impossible  to  construct  a truly 
interelement  compatible  finite  element  from  this  set  of  solution  functions 
without  resorting  to  subdomain  techniques  (Tocher,  Ue  Veubeke,  ref's  6,7) 
Moreover,  obtaining  interelement  compatibility  in  Pian's  sense  (ref  1 ) 

is  not  dependent  upon  polynomial  representation  in  the  interior. 
Furthermore,  there  exists  no  recursion  formula  that  allows  finding  in 
an  expedient  manner  successive  terms  of  the  series  of  independent 
polynomial  solution  functions  to  the  homogeneous  biharmonic  equation. 
Therefore,  such  a series  is  not  particularly  well  suited  for  the  develop- 
ment of  a family  of  finite  elements  with  an  arbitrary  number  of  nodal  points. 
Furthermore,  it  is  generally  accepted  that  polynomial  series  representations 

of  general  continuous  functions  in  the  majority  of  cases  show  poor  con- 
vergence when  compared  with  trigonometric  representations  of  the  same 
functions.  In  other  words,  a smaller  number  of  parameters  is  needed  for  a 
given  maximum  error  when  a smooth  function  is  represented  by  a tri- 
gonometric representation  as  opposed  to  a polynomial  representation. 

Since  it  is  required  to  find  the  inverse  of  the  matrix  S (eq.4)jit  is 
desirable  to  keep  the  size  of  this  array  as  small  as  possible  for  a 
given  refinement  of  the  finite  element.  This  speaks  clearly  in  favor  of 
trigonometric  representation  of  the  displacement  function.  Any  savings 
in  computer  time  in  the  inversion  process  will,  however,  be  considerably 


1/1 


degrees  of  freedom 
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leduced  since  the  evaluation  of  each  trigonometric  and  exponential  term 
performed  on  a digital  computer  itself  involves  evaluation  of  truncated 
polynomial  series.  Foregoing  remarks  indicate  that  the  best  choice  of 
the  functional  representation  of  the  transverse  displacement  of  the  plate 
for  the  equilibrium  finite  element  with  a large  number  of  degrees  of 
freedom  is  not  obvious.  The  set  of  solution  functions  to  the  biharmonic 
equation  that  is  chosen  here  are  the  Le^  series  solutions  for  the  refer- 
ence field  with  various  boundary  conditions.  Donaldson  et  al  (refs. 8 thru  12) 
superimpose  the  Navier  solution  and  Le\^  solutions  for  the  dynamic 
plate  bending  equation  for  the  purpose  of  finding  the  response  to 
harmonic  excitation  and  for  finding  the  eigenvalues  for  the  rectangular 
trapezoidal,  and  triangular  plates  for  various  boundary  conditions. 

This  method,  the  Extended  Field  Method,  embeds  the  plate  to  be 
analyzed  (the  actual  plate)  in  a larger  rectangular  plate  (the  extended 
plate).  The  Lei^  solutions  1,4  , chosen  are  those  that  satisfy 

simply  supported  boundary  conditions  on  all  but  the  i^^  edge.  The 
arbitrary  transverse  displacements  and  rotations  of  the  i^^  edge  of 
W^{%)  are  determined  by  the  free  parameters  of  the  series. 

The  boundary  conditions  of  the  actual  plate  are  then  enforced 
in  an  approximate  sense  by  requiring  the  error  to  be  orthogonal  to 
a set  of  test  functions.  This  is  precisely  the  technique  proposed 
here,  (eq.4),  for  finding  the  coordinate  transformation  for  a displacement 
compatible  finite  element. 

However,  the  Navier  and  Levy  series  solutions  to  be  used 
to  produce  a finite  element  stiffness  matrix  and  the  associated 
force  vector  are  those  that  satisfy  the  static  equation.  The 
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appropriate  Navier  solution  is  given  in  ref. 5 and  below.  The 
appropriate  Levy  series  solutions  are  derived  in  the  sequel.  It 
nay  be  noted  here,  that  it  is  customary  in  finite  element  solutions 
of  dynamic  problems  to  use  mass  matrices  that  are  based  on  a smaller 
nunber  of  inertia  references  than  displacement  coordinates  used  for 
developing  the  structural  stiffness  matrix.  The  stiffness  matrix  is 
then  condensed  and  the  solution  of  the  dynamic  problem  is  found  from 
the  reduced  system  of  equations.  This  is  justified  by  the  fact  that 
the  solution  of  the  dynamic  problem  is  more  sensitive  to  the  refine- 
ment of  its  elastic  references  than  to  the  corresponding  refinement 
of  its  inertial  references. 

In  order  to  produce  a finite  element  stiffness  matrix  it  is 
necessary  that  the  displacement  function  used  contains  components 
that  permit  elastic  comer  displacements  and  rigid  body  displacements.  This 
can  be  achieved  by  the  use  of  a few  polynomial  terms.  The  completeness 
requirement  for  the  plate  bending  element  requires  for  the  set  of  poly- 
nomial terms  that  at  least  all  terms  up  to  the  second  order  are  present 
in  the  displacement  function.  This  means  that  in  addition  to  all  rigid 
body  displacements,  all  constant  strains  should  be  included.  The  elastic 
comer  displacements  can  be  represented  by  the  polynomial  function 

W.M)-C,(l-S)t  » C.(l-I)(l-J)  1-C,  . (13 

It  is  observed  that  in  addition,  these  functions  are  able  to  represent 
all  rigid  body  displacements  and  constant  twist.  Unless  the  number  of 
LeA^  series  terms  employed  is  sufficiently  large  to  approximate  constant 
bending  strain  well  enough,  it  will  be  necessary  to  include  in  Wo(^) 
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the  functions 


The  solution  for  the  inhomogeneous  biharmonic  equation 

V‘‘W 

subjected  to  simply  supported  boundary  conditions  at  the  edges 


X-  0 , 


a 


, ro  . 


of  the  reference  field  (the  extended  plate)  is  given  by 

p-.sin'f  sinf*. 

fn®l  n»i  ’ 

where  . b 

d%  . ,14) 

The  Levy-series  solution  for  the  extended  plate  with  a simply  supported  boundar>^ 
along  all  edges  except  for  example  the  edge  X»Q  can  be  found  by  choosing  for 
the  displacement  function  the  Levy  series 

Wj()»0  - l' X,(J:)sinfJ  . 

n*i 

The  index  3 indentifies  the  edge  of  the  extended  field  at  which  arbitrary 
displacements  and  rotations  are  prescribed  (Fig.  2 ).  This  series  satisfies 

the  boundary  conditions  at  the  edges  ^^0  and  automatically.  Substitu- 

tion of  this  displacement  function  into  the  homogeneous  differential  equation 
and  observation  of  the  orthogonality  ofl/W^over  the  interval  (O^b)  yields  the 
ordinary  differential  equation 

x:  -2fa:  - o 

whose  solution  is  subject  to  the  boundary  conditions  XM-Xjol’  0, 
X.(a)-h,  and  X„  (a)-  k«.  Substitution  of  the  exponential  function  6 
into  the  equation  and  cancellation  of  the  ccmmon  factor  gives  then  the  char- 
acteristic equation  T - “ 0 

which  has  the  double  roots  f"!  The  solution  of  this  differential 
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equation  is  therefore 

X.(«f)  - Ae'"’  + 5,6'*"'  -I-  r(Cs^'  + Se*"*) 

or  equivalently 

X,(x)“A„co5/iA.x  +Ei,smliA„x+CA^“5^>^-^+  D.^xsinhA,x  ^ 


. A = 

where  /\n  u • 


Aiy-:liCu1:ion  of  the  boundary  conditions  yields 


r 

1 

0 

0 

0 ][/ 

u 

0 

1 

0 

0 

Z E 

• 

0 

cosh  ia 

sinhio 

Aacoshio 

AiflsnhAna  C 

1 

h. 

.sinhifl 

co5hAo 

[oshio+iosMnCi  sinh/\na-^^acoshiQ|{[ 

1;  1 

t 

Solving  this  matrix  equation,  the  coefficients  are 

An  “ Dn  “ 0 

Bn  - 3.(|)  ( -hn^shiQ  f AnQsinfilo)  iacoshla  } 

and 

C,  - ].(§)i  Lcoshia  -^sinhAnfl  j , 


where 


1(f)-  [A. a - sinh Anfl cosh Anfl ] ' . 

Substitution  into  the  expression  forW^,  eq.  ( 15  ), gives 

WjM  “ I ’j,j(^)5inA..^{h,[-(coshX,a+A,.a5inhA«a)s*nh?LJ:+cosh/\,Q  -/IxcoshA-r  J 

+ [ A,a  cosh  id sinh  U - sinh  ka  kx  cosh  kx]] 

Am 

Similar  solutions  can  be  found  for  W,,  » i “ 1,2  and  4. 

At  this  instant  it  should  be  noted  that  ) is  a constant  for 

given  n and  ordered  side  ratio.  Since  there  is  no  need  to  attach  a specific 
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meaning  to  the  parameters  hn  and  , convenient  scaling  factors  may 

be  chosen  that  simplify  the  algorithm  to  be  developed  and  avoid  difficulties 

a 

that  might  arise  from  computing  large  numbers  of  the  form  C 

In  previous  work  on  the  EXFM  ^vstiich  dealt  strictly  with  the  dynamic 
plate  equation  the  corresponding  factors  Jui^)  had  singularities  for  certain 
side  ratios  ^ and  n . Therefore  such  simplifications  \jvre  not  possible. 
Thus  let, 


Icosh  AjA  • ^coih  lijX  - [ cosh  A, A + sinfi  A, A ]5inh  AjYl 

and 

Kj(X.Y;A,B)  = e-^"sinA,r  - ,i6) 

{ - sinh  A/  • A,X  (o:h  A,X  + X,A  cosh  A, A sihhAjA } , 

where 

Aj  - jJT/B 

With  this  definition,  the  displacanent  functions  that  correspond  to  simply 

supported  boundary  conditions  along  all  edges  of  the  extended  plate,  but  the 
i^^  edge,  where  arbitrary  displacements  are  prescribed,  become 

~ I h.j  Kj  (a-x^;  Q,b  ) “I  kij  Kj(a-x,i^  ;a,b) 

Wi(X)  « I h^i  lfj(  y,x;  b,Q  ) + Z * k^XjC  ; b,a  ) (i?) 

wjx)  - I'  b,a)  - J]  . 

The  transverse  displacement  of  the  plate  can  thus  be  expressed  by 

* Wo(X)  -b  Z W//X)  +WpU).  1.18) 

To  obtain  the  displacement  of  the  plate  edge  in  terms  of  the  displacement 
parameters  for  the  plate  interior  it  is  simply  required  to  evaluate  the 
displacement  at  the  edge  points.  Using  edgewise  boundary  coordinates,  the 
transformation  from  boundary  coordinates  to  reference  coordinates  is  given 
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by  , IJ". , 

where  are  the  ccxjrdinates  of  the  comer  nodal  point  i and  the  ratios 

. ajid  ^ give  the  directions  of  the  straight  edge  boundary  segment  i. 
For  representation  of  the  total  boundary  by  one  continuous  parameter  this 

transformation  becomes 

where  - 0 

To  obtain  the  edge  rotation,  it  is  necessary  to  first  take  the  derivative  in 
the  direction  of  the  outside  normal  and  then  evaluate  this  derivative  at  the 
edge  points. 

The  transverse  displacement  and  rotation  at  the  plate  edge  thus  expressed 
in  terms  of  the  parameters  of  the  interior  representation,  it  remains  to 
express  these  displacement  quantities  in  terms  of  the  nodal  degrees  of  free- 
dom. The  transverse  displacement  and  rotation  of  the  finite  element  bound- 
ary have  been  chosen  (sect.  2.1)  to  be  splines  in  terms  of  the  nodal  degrees 
of  freedom.  For  the  edge  displacements  of  edge  j the  cubic  spline  is 

where 

!,(«)-  I 1 

] w.., 

with  A,  - 

and  »'•!  j and  uefincri  by  figure  3. 

Here  Id)  is  the  Heavyside-Unit  step  function  and  the  edge  coordinate 
at  nodal  point  i , ITlj  the  nodal  point  number  corresponding  to  the  comer  node 
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of  corner  j . 

The  linear  spline  for  the  edge  rotation  at  edge  j is  given  by 

■ tm'  1 !(«-«,)-  , 

where 

9^(e)  - [ 

with 


r20) 


2.3  The  Rectangular  Plate  Bending  Finite  Element 

The  experiences  gained  in  references  (9,10,11)  can  be  applied  to  the  pre- 
sent task  of  developing  a rectangular  finite  element  for  plate  bending  via 
tnis  equilibrium  approach.  First  of  all,  there  is  little  reason  to  have  the 
actual  rectangular  finite  element  smaller  than  the  rectangular  extended  plate. 
Secondly  computational  advantages  are  gained  if  the  error  of  the  displacement 
and  rotation  is  required  to  be  orthogonal  to  a set  of  test  functions  for  each 
edge  separately.  That  is, ^instead  of  equation  ( 4 ) the  boundary  error 

reduction  is  given  by  | Grr0r'(^)<^(^)d^  = 0 , = j . The  total 

number  of  such  equations  must  again  equal  the  number  of  parameters  in  the 
displacement  functions  for  the  plate  so  that  the  coefficients  matrix  associated 


with  these  parameters  is  square  and  nonsingular.  Also  some  simplification 
in  the  equations  for  the  force  quantities  are  noted.  It  would  be  quite  simple 
to  derive  the  Force-matching  conditions  directly  in  terms  of  derivatives  in 
the  X and  directions,  the  reference  coordinates,  without  reference  to  edge 
coordinates  which  require  transformation  coefficients.  However,  to  maintain 
clarity  and  uniformity  the  structure  of  the  development  is  retained  and  the 
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simplification  of  the  transformations  are  merely  noted.  The  transformation 
coefficients  , y^aJk  and  vanish  identically  and  the  remaining 

coefficients  are  given  in  table  1 


Transformation  coefficients  for  the  edge  coordinates  for  the  rec- 

tangular plate. 


\f? 

N 

1 

2 

5 

4 

2 

D(i-p)(-2) 

D(I-p)(2) 

D{l-P)(-2} 

D(i-p)(2) 

5 

D 

0 

-D 

0 

0 

D(2-p) 

0 

5 

D(2-p)(-i) 

0 

D{2-P) 

0 

6 

0 

-D 

0 

D 

7 

-D 

0 

D 

0 

9 

-Di; 

0 

Du 

0 

10 

0 

-Pp 

0 

Dp 

12 

i 0 

-D 

0 

D 

Table  1 

24 


2.4  The  Quadrilateral  Plate  Bending  Finite  Element 

To  obtain  the  transformation  matrix,  equation  5,  from  the  field 
displacement  coordinates  to  the  nodal  degrees  of  freedom  one  may  require 
the  displacement  discrepancy  to  be  orthogonal  to  a set  of  test  functions 
for  each  edge  separately  or  for  the  closed  boundary.  In  either  case  the 
total  number  of  resulting  equations  must  equal  the  number  cf  tlie  [Brameters  of 
the  displacement  function.  The  alternative  procedure  preferred  by  Donald- 
son and  Chander,  (refs.  8,9),  in  which  the  intervals  over  which  the  error 
is  made  orthogonal  to  the  set  of  test  functions,  extend  to  the  edges  of 
the  extended  plate  is  certainly  attractive  for  the  case  of  (determined) 
homogeneous  boundary  conditions  since  computational  advantages  arise. 

For  the  purpose  of  developing  a stiffness  matrix  for  a finite  element 
where  the  boundary  conditions  on  the  finite  element  remain  undetermined, 
extension  of  the  intervals  for  the  orthogonalization  process  would  add 
arbitrariness  to  the  process  and  with  it  an  additional  source  for  error. 

Chander  also  observed  deterioration  of  convergence  as  the  extended 
field  becomes  larger  in  comparison  to  the  actual  field.  One  should  there- 
fore choose  the  extended  rectangular  plate  such  that  the  area  of  the  con- 
vex quadrilateral  plate  is  as  near  as  possible  in  size  to  the  area  of  the 
extended  plate. 


2.4a  The  Coordinate  Transformation  for  the  Quadrilateral  Plate  Bending 
Finite  Element 

The  transformation  (eq.  5)  that  relates  the  nodal  degrees  of  freedom, 
the  , to  the  parameters  hj  of  the  displacement  functions  for  the  plate 


interior  results  from  matrix  operations  on  equation  (4).  This  latter 
equation  is  obtained  by  demanding  that  the  discrepancy  in  the  external 
and  internal  representation  of  the  same  edge  displacement  quantitites 
be  orthogonal  to  some  set  of  test  functions.  Care  must  be  taken  in 
selecting  this  set  of  test  functions  so  that  the  matrix  0)  in  equation  (5) 
is  square  and  invertable.  For  each  of  the  four  edges  the  discrepancies  in 
the  transverse  displacement  and  in  the  rotation  are  both  required  to  be 
orthogonal  to  an  orthogonal  set  of  functions.  A convenient  orthogonal  set 

vlTg 


is  given  by  for  each  edge 

simultaneous  equations  is  given  by 


The  resultant  set  of 


ftk 


t sin^  dg;  »0  j-l,N2k  ;k-l,4 


(21) 


Substitution  of  the  displacement  functions  into  these  equations  gives  the 
four  sets  for  the  edge  displacement 

A 


+ sinf^de 

A 

c,  1 [i  ( de, 

It, 

+ C»  [[^(7n,,-h  ^)I [ i ( ^«k  ^ 1 

+ Cs  [(a)  (4k  + 2xmk(^),^-b(^)V  ) Sinf^d^ 
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27 


where  A;  , 


and  the  four  sets  for  edge  rotation 

+C,  l(^)h  Sin^fd? 


t 
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+ de;  } 

^"iiTr 

1-1...  y l<”  Ij^ 


(22b) 


zx  means  differentiation  of  the  functions  KandX  respectively  with 
respect  to  the  first  argument,  while  r.x  means  differentiation  of  the 
functions  ]{  andV  respectively  with  respect  to  the  second  argument. 


t 


2.4b  Virtual  Work  and  the  Stress  Resultant -Nodal  Force  Relation 


The  transformation  from  the  nodal  degrees  of  freedom  to  the  para- 
meters of  the  displacement  functions  for  the  plate  interior  is  given  by 
the  simultaneous  set  of  equations,  equations  (22a)  and  (22b)*  of  the 
previous  section.  It  remains  to  obtain  the  elements  of  the  stiffness 
kernel  and  the  coefficient  matrix  of  the  loading  parameters  in  the 
expression  for  virtual  work.  The  virtual  work  expression  (eq.lOjcan  be 
written  in  the  form 

^ k=i  « • > (23) 

^)\  d^q, 

where  ]fp  are  here  (generically)  the  displacement  functions,  i.e.  the 
polynomials  representing  rigid  body  displacements,  elastic  comer  displacements 
and  constant  strains  as  well  as  the  hyperbolic- trigonometric  functions  that 
represent  the  remaining  elastic  solution  functions  of  the  homogeneous  equation. 
The  ^ are  the  displacement  functions  corresponding  to  the  inhomogeneous 

7,  ‘i’ 

equation,  the  CU  represent  the  appropriate  partial  differential  operators 


in  equation (23) 


Jil  ~ i , 1-1/ 

Z)i  - ^ j 1-I0J2 


, ^ - 7,9 


rl"  /x- 

, (^2  * wSij 

, , A’-  ^ 


iZ 


and  the  vectors, hi  , (I|  p are  given  by 


f ' 

c. 

w. 

'p,. 

u 

Pu 

h.) 

piNp 

hm. 

Jl  = ■ 

■ 

P = ‘ 

pll 

k.. 

kiK. 

b 

‘ PmrNp, 

To  economize  an  abbreviated  notation  is  introduced  that  reduces  the  poly- 
nomial solution  functions  to  a single  form  and  the  hyperbolic  solution 
functions  to  a single  form  with  four  different  sets  of  arguments.  This 
gives  for  the  polynomial  terms 

W.(X)  -i  a;  + a'*  + (24) 

where  the  appropriate  coefficients  of  the  polynomial  terms  are  collected 
in  table  2 . 

The  coefficients  Qj  for  the  polynomial  solution  functions 


I 

z 

3 

4 

5 

6 


0 


0 

I 

0 

0 

0 

0 


0 

a 


1 

a 


2 3 4 5 


b 0 sis  0 

■ b 0 ab  0 

0 0 ~^b  0 

0 0 5B  0 

0 0 0 
0 0 0 


0 

0 

0 


Table  2 


For  the  hyperbolic- trigonometric  functions  , A ,B  >*nd  A,B) 

• 3 

the  new  notation  is 

cfy  (x,Y,  r,s)  . 

where  the  integer  parameters  r and  s are  to  be  interpreted  as  follows: 

^ r - 1 

and  * 3K  r * Z 

and  the  integer  s matches  the  first  subscript  of  the  distributed  parameters 
h.  and  kii  with  the  additional  meaning  that  if  r = 2 and  s = 1 or 
s = 4 the  negative  of  the  function^  should  be  taken. 

Thus, 

(X,Y;A,&;05}  “ AjY(i4  AjXC0sh?lj>(+ B Sinh  AjX)/25] 

where  Aj  “jH/B 

and  ^ cosh A , S = "[cosh^A-t-AjA sinhAjA ] r«l,s*i,4 

A-  « sinh  AjA  ; B • " ^jA  cosh 7\jA  when  r-z 

= - smh  AjA  , S*A^Aco5hA^A  • hen  r=z  s-{2,3} . 

It  is  convenient  to  define 

bl*"  - a:  + a:  x,.,  + a\  + q\x\  + a\  t a\ 

bf-  a:(^p:  + a:(dp^id^\  + 
til:  - ai  +zaU-.  + ajij,. 
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b',.‘  - 2a:(^).  + 


bl,  ’ a:  + a;x^  + 2a‘,i|., 
b'y*  a'.^l  +2a'.(^), 


so  that 


K >S' 


\No yVo ,«(^)“Z 2 Os  C f 
'•  >i*' 

4 

W»  iVo  O'*  C X 

Wo  (g,]  =lcb:'+  + b:^^’')c, 

^■1  ' 

Wo.  bix^)  C/,  ^ 

-'  X'' 

6 i k i k 

and  “Z  (boy  b„j^)C,  ^ 

Using  these  definitions  a few  typical  coefficients  of  the  matrices 

llj"  ‘[l)j}':Mf)(cdTX^)  dsf 

““  19J  -[l\iflx,){fft)ckf 

are  worked  out  and  a general  scheme  is  provided  in  the  appendix. 

Of  course,  the  stiffness  kernel  and  the  matrix  IR  in  equation  (11)  are 


given  by 


xal  k«l 


R=il 


^1  k<l 


5.  Appendix--  Matrix  coefficients  for  the  Stress  Resultant  - Nodal  Force 
Relation.  ^ 

Components  of  the  matrices  [W 
The  typical  components  are  of  the  form 


I,  12  ; k-  1^11  ^ 

/■i'i' 

where  the  are  the  differential  operators  defined  in  section  2. 4b, the 

d^n  and  are  the  solution  functions  for  the  complimentary  and  particular 

solution  of  the  differential  equation  respectively. The  index  H in 

I) 

ranges  from  1 to  6 + Z (Hj+Kjj.  In  the  subranges  of  H : {l,  ,6}  ^ 

f7.  ,6+H,l  , {7+Hl,',6'+H.+K,i,  ..  (7-|(H,kK,)+H.,  ,6+^_(H(+K,} 

the  function  Hn  takes  different  forms. The  definition  of  for  the 
different  subranges  is  given  below  in  table  A1  . The  required  derivatives 
of  the  polynomial  functions  are  given  in  section  2.4b  , those  for  the  hyper- 
bolic-trigonometric function  ^ are 


/’•  A,  sinAjY  [£DshA,X-fA,Xsinl)AjJ(l  + SjCoshA/ t 
r • \ msAjY  M A,AcodiX,X  * B,sinhW  1 

a;  OTA,Yl41?5inhX.X-A,XcosdAXl +i,,an  A/l 

r.  a;  e-^*‘''a)sA,YMlMdiykA,Xsinhyl  +SjC»hAj)(  i 
sitiA,Yf4  A,Xc()diX,X  sinM,X  1 

• Aj  sinA,Y  {^[3aishAjX+A,XsirihX/l  + Sj  os/iXjX  i 

/”'•  A;  osAjY  f-^IXsinhytAi/cashAiXl +B|  anhA/ i 

r"--A’  sinAY  f«^kos)iAA  + AXsitiAj(  1 * I 

cosA,Y  i4  A.XoiJhAjX  tS,siiHX( 
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Range  of  n 


Fund  I on 


Range  of  index  j 


I to  6 

7 to  6n-H.  ^(Q-x^ij- a,b;  Ij  I ) 

7+  H.  to  6+  H,+K.  ifj  Q,  b ; 2^  I ; 
r+H.+IC,  to  i^^z-ha;l,2) 

7+iH,+K  to  fj(  y^x;b^Q;7^2^ 

7+lW^^.)  to  k :(,u;Q,b;2,3) 

to  6+J(H.tfCj  <:§(  x,t^;a,b;2,5) 
A^Oi+ki)  to  6+?(H,^■l(i>■H^  ^(b-i^/;b,Q; 

7+|ffi+KM  to  6t|(H,fl(4)  Jj(b-ij,x)b,Q}2A) 


]= 

-.6 

j* 

. ,H 

J- 

i, 

r 

I,. 

3= 

V 

i’" 

i>- 

r 

i,. 

,K, 

j' 

I, 

i- 

r,s)  * e ^ sin^jY(s^  co5fi/^i](  •+^5inhAjX) 
where  = jlT/S 

and  5z^=  Cosh  Aj4  , ^^  = - [ Cosh  AjA  +AjA5inhAjAl^w;hen  / ■ I ^ s =1;4 
sinh  Aj  A , - AjA  cosh  AjA  r*  2 , s^l  i,4} 

4=*sinh^A,  0^=  ^Acosh!^jA  r«z,v[2,3} 

ihe  cons1$nfs  q"  are  given  in  table  2 of  section  2.4  b 

Table  A1 

Solution  Functions 


With  these  definitions  and  those  for  the  coefficients  in  section  2."b 
it  is  a simple  matter  to  vrrite  out  the  coefficients of  the  matrix  IH 
Example  1 

The  elements  of  the  matrix  IH  derive  from  the  term 

[ W,xx(X«,)  - 

of  equation  ( 10  sect.  2.1). 

Thus  for  p ^ \ f A 

Jfpq”  2 (al  -a’Jbl'' 

in  particular  for  p **  3 ^ Q - 2 


4k 


3fv- 

0- 

1 4 k =• 

b. 

0 

or  for 

p ^3 

q-5 

1 

(a'f”. 

_ -L 
‘ Ob 

for 

p-  z 

6+H. 

q= 

5 

t.- 

2 

5* 

1 

;QA;  I;)) 

^ e AjQ  • Ai(a-xj'Cosh 

-[coshA,a+Aiasir]hA,alsinhAj(a-;c«,) } ,j= 


with 


A - -f  ^ 
J b 


Example  2 

The  elements  of  the  matrix  H derive  from  the  term 

of  equation  (10  sect.  2.1  ) 
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f 


For  p ,C^  ■ I , ,6 


3fp,  = 0 

For  p * 1, , 6 and  7+H.,,6+Hi‘*'K, 

Jfpf  2,i)(bI+b''^+b';V)cl^  ^-i  = b,K. 

= - fsin  sinhAjQ  [3coshAj(a-T.:^{4)+A,lQ-r4^\? ) 

siiihAj(a-r.;(^^e)“AjQcoshAia  5inhA,(a-x.,-^^i^]l  d§ 

, nrr^  ^ 


where 


specializing  further  to  P''6  then  Ol^'O  except  <3=^  so  that 

K + b!^+b‘ V = • 

Example  5 

The  elements  of  the  matrix  IH  derive  from  the  term 

L 

' 

of  equation  (10  sect.  2.1  ) 

For  P Z ,6+H,  and  p = 7 +I  (H/K,)  , , 6 + 2(H.  + KJ  + H., 

^ <•!  >■! 

IP  = -2(;\»btA.Q)  ( ^ ^ , n , 


1 pf  Xm  ^r  ^ I Sj  b,a ; i, 4 ) 

■ i'ca-PCviti? ; ; Q,t3;  >,  1 ) de^  ^ 

\ mTT  * Jill' 

where  Aw  ” and  A,  (3 

W - \,  , H4  and  ; H,  , 

The  coefficients  of  the  matrices  are  obtained  in  similar  manner,  here 
the  index  p runs  from  I to  Mo^  No .The  function  ^ in  equation  (25  ) has 


the  form 


f^i%)  - Sin  Sin  ^ ^ 


L 


where  following  correspondence  between  the  index  CJ  and  the  index  pair 
(m,n)  is  set  up 


i 


?; 

V- 

h 


q m . n 

1 > ; I 

2 z , 1 

M.  M.  , I 

M.^1  I , z 

k M,  , 2 

2Mo+l  I , 3 

M.xN.  M.  , N. 


Exanrple  4 

The  elements  of  the  matrix 


derive 


as  in  example  2 


from  the  term 


A 

de, 

For  p « 6 and  q = 5Mo+3 


HO 


4.  Appendix--  Summery  of  the  Extended  Field  Method  of  Analysis 
as  applied  to  Uniform  Thin  Plates. 

The  governing  differential  equation  for  the  harmonic  response  of  thin 
uniform  plates  is 


D VV  - yc*3*W  » piX)  Sl’n  Oit  ^ 

where  WCXj't ) ■=  W(X)Sinaj^  is  the  time  dependent  transverse  displacement 
of  the  plate,  D is  the  bending  rigidity,  CO  is  the  excitation  frequency 
and  pCX)  is  the  spatial  distribution  of  the  transverse  pressure  , and  ^ 
is  the  mass  per  unit  area.  The  amplitude  1A/(X)  of  the  transverse  displace- 
ment can  be  expressed  as  the  sum  of  the  solution  to  the  inhomogeneous 
biharmonic  equation 

Dv'W  - Y<a*W  = pix) 

for  a rectangular  reference  field  (the  extended  plate)  with  all  edges 
simply  supported  (the  Navier-solution)  which  is 

and  the  four  Levy-series  solutions  , for  the  same  extended 

plate  with  simply  supported  boundary  conditions  on  all  but  the  i^^  edge 
''’hich  undergoes  arbitrary  displacements  that  are  characterized  by  the  para- 
meters h;j  and  kij  .For  example  the  Levy-series  solution  for  the  1'"'^  edge 
is  given  by 

W,(x.)  = ? J I h.j[siCOS(h)s,Q  iinhij(a-a:)  - rjCoslir^a  iiii(]i)S)(a-x)  1 

+ kj  [sin{}|)SjQSinhr^(a-x)  - 5inhr<QSin(lT)Sj(Q-z)]  i , 

where 

3^-  [Sj  sinh  r^Q  cos(h)s^Q  - r^ooshfjQ  sinchj^a  J~‘ 
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with 

r,  - [ if)’ + 4)(^)'‘l'' 

5|  - I 1*  1<  i 

Here 

r 

and 


Sin(h)s,  jc  *' 

sinh  SjX 

cos(h}  Sjjr  = ' 

cosh  ^x 

The  total  displacement  of  the  plate  can  then  be  written 

<5 

W(x,4)  = sin  c5l)4  • Z W^-(x) 

The  solution  to  the  harmonic  excitation  problem  of  an  arbitrary  quadri- 
lateral (convex  region  ) plate  that  lies  within  the  extended  plate  with 
any  combination  of  clamped  and  simply  supported  edges  and  perhaps  a single 
free  edge  is  then  obtained  by  approximately  enforcing  the  desired  boundan' 
conditions  at  the  edges  of  this  embedded  actual  plate.  In  the  case  of  a 
free  comer  , i.e.  two  adjacent  free  edges  , the  solution  function  must 
be  augmented  by  a comer  solution  which  , for  example  , for  the  edges  2 
and  4 being  free  is  given  by 

- ^ 

Wo(x)  **  c I sinh(j>(Q-it)sir)h«f(b-^)+sinf(a-x)sin(^b-i^)l  , 

The  edge  numbering  corresponds  to  that  of  figure  1 of  section  2.1 


i 


The  solution  to  more  complex  structures  such  as  beam  stringer  supported 
plate  systems  can  be  obtained  by  approximately  enforcing  the  matching 
conditions  along  common  boundaries  (ref.  8 ).  The  set  of  simultaneous  equa- 
tions that  determine  the  parameters  and  kij  is  obtained  by  requiring  the 
error  at  the  boundary  edge  or  matching  edge  to  be  orthogonal  to  some  set  of 
independent  test  functions.  The  EXFM  has  also  been  successfully  employed  in 
finding  the  eigenfrequencies  for  uniform  quadrilateral  plates  with  various 
boundary  conditions.  Chander  employed  the  method  of  false  positioning  to 
extract  the  eigenvalues  which  appear  in  a transcendental  manner  in  the 
matrix  equations  that  result  from  the  process  of  error  orthogonal izat ion. 
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